Write a programme following the method of Monte Carlo Integration to calculate
$$ I = \int_0^2 \sin^2 [\frac{1}{x(2-x)}] dx. $$As you will need to calculate $f(x) = \sin^2 [\frac{1}{x(2-x)}]$ many times please write a user defined function for this part of your programme.
In [6]:
import numpy as np
import matplotlib.pyplot as plt
#This just needed for the Notebook to show plots inline.
%matplotlib inline
#Define the function
def f(x):
fx = (np.sin(1/(x*(2-x))))**2
return fx
#Integrate the function from x=0-2
#Note that you need to know the maximum value of the function
#over this range (which is y=1), and therefore the area of the box
#from which we draw random number is A=2.
N=100000
k=0
for i in range(N):
x=2*np.random.random()
y=np.random.random()
if y<f(x):
k+=1
A=2.
I=A*k/N
print("The integral is equal to I = ",I)
In [7]:
#Calculate the error:
sigmaI = np.sqrt(I*(A-I))/np.sqrt(N)
print("The integral is equal to I = ",I)
print("The error on the integral is equal to sigmaI = {0:.4f}".format(sigmaI))
In [8]:
#Draw N values from the distribution, and calculate their mean to
#use the mean method for integration.
xvalues=[]
fvalues=[]
for i in range(N):
x=2*np.random.random()
y=np.random.random()
if y<f(x):
fvalues.append(f(x))
xvalues.append(x)
fmean=np.mean(fvalues)
Imean = 2*fmean
Imeansigma = 2*np.sqrt(np.var(fvalues))/np.sqrt(len(fvalues))
print("Mean Method: ")
print("The integral is equal to I = {0:.4f}".format(Imean))
print("The error on the integral is equal to sigmaI = {0:.4f}".format(Imeansigma))
print("The percent error is: {0:.2f} ".format(100*Imeansigma/Imean))
print("**********************")
print("Compare to the `hit or miss` Monte Carlo Method: ")
print("The integral is equal to I = {0:.4f}".format(I))
print("The error on the integral is equal to sigmaI = {0:.4f}".format(sigmaI))
print("The percent error is: {0:.2f} ".format(100*sigmaI/I))
In [9]:
#Using mpl_style file.
import mpl_style
plt.style.use(mpl_style.style1)
plt.plot(xvalues,fvalues,'.')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.show()
In [18]:
np.mean?
In [ ]: